Load libraries

if (!require(R.utils)) install.packages('R.utils')
## Warning: package 'R.utils' was built under R version 4.0.5
library(R.utils)

if (!require(dplyr)) install.packages('dplyr')
## Warning: package 'dplyr' was built under R version 4.0.5
library(dplyr)

if (!require(stringr)) install.packages('stringr')
library(stringr)

if (!require(Seurat)) install.packages('Seurat')
library(Seurat)

if (!require(purrr)) install.packages('purrr')
library(purrr)

if (!require(scCATCH)) install.packages('scCATCH')
library(scCATCH)

if (!require(ggplot2)) install.packages('ggplot2')
## Warning: package 'ggplot2' was built under R version 4.0.5
library(ggplot2)

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!require(SingleR)) install.packages('SingleR')
## Warning: package 'matrixStats' was built under R version 4.0.5
## Warning: package 'BiocGenerics' was built under R version 4.0.5
## Warning: package 'GenomeInfoDb' was built under R version 4.0.5
library(SingleR)

if (!require(msigdbr)) install.packages('msigdbr')
## Warning: package 'msigdbr' was built under R version 4.0.5
library(msigdbr)

if (!require(tibble)) install.packages('tibble')
## Warning: package 'tibble' was built under R version 4.0.5
library(tibble)

if (!require(fgsea)) install.packages('fgsea')
library(fgsea)

#if (!require(clusterProfiler)) install.packages('clusterProfiler')
library(clusterProfiler)

#if (!require(org.Hs.eg.db)) install.packages('org.Hs.eg.db')
library(org.Hs.eg.db)

if (!require("org.Mm.eg.db", quietly = TRUE))
    BiocManager::install("org.Mm.eg.db")
library(org.Mm.eg.db)

Obtain files

Download data

This project uses the data from Grosselin et al, 2019, which is uploaded to GEO.

mkdir -p data
cd data
wget -r -np -nd 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE117nnn/GSE117309/suppl/GSE117309_RAW.tar' -R "index.html*" # Omitting the index file

Decompress files

# Untar the main file
untar("data/GSE117309_RAW.tar", exdir="data/")

# Untar sc-RNAseq files
RNA_files <- list.files(path="data/",pattern="scRNA")
untar_RNA_files <- function(file){untar(paste0("data/",file), exdir="data/RNA_files")}
sapply(RNA_files, untar_RNA_files)

Read files into R

# Create function to load the files
load_scRNA_files <- function(model, organism){
  directory <- paste0("data/RNA_files/filtered_gene_bc_matrices_",model,"/",organism,"/")
  expression_matrix <- ReadMtx(
  mtx=paste0(directory, "matrix.mtx"),
  features=paste0(directory, "genes.tsv"),
  cells=paste0(directory, "barcodes.tsv"))
    # Correct gene names for easier automatization
  expression_matrix@Dimnames[[1]] <- str_remove_all(expression_matrix@Dimnames[[1]],
                                             pattern=paste0(organism,"_"))
  seurat_object <- CreateSeuratObject(counts = expression_matrix, project=model)

  return(seurat_object)
}

models <- c("HBCx-95","HBCx-95_CAPAR", "HBCx-22","HBCx22-TAMR")

scRNA_files_hg19 <- mapply(FUN=load_scRNA_files, models, organism="hg19")
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
names(scRNA_files_hg19) <- paste0(names(scRNA_files_hg19), "_hg19")

scRNA_files_mm10 <- mapply(FUN=load_scRNA_files, models, organism="mm10")
names(scRNA_files_mm10) <- paste0(names(scRNA_files_mm10), "_mm10")

Analysis of human cells

Data integration

# Normalize and identify variable features
scRNA_files_hg19 <- lapply(X = scRNA_files_hg19, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# Select features repeatedly variable between the 4 datasets
features <- SelectIntegrationFeatures(object.list = scRNA_files_hg19)

# Select anchors for integration
brca_anchors <- FindIntegrationAnchors(object.list = scRNA_files_hg19, anchor.features = features)
## Warning in CheckDuplicateCellNames(object.list = object.list): Some cell names
## are duplicated across objects provided. Renaming to enforce unique cell names.
## Scaling features for provided objects
## Finding all pairwise anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2276 anchors
## Filtering anchors
##  Retained 1594 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2558 anchors
## Filtering anchors
##  Retained 1046 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 3666 anchors
## Filtering anchors
##  Retained 1570 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2201 anchors
## Filtering anchors
##  Retained 1215 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2766 anchors
## Filtering anchors
##  Retained 1635 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 3104 anchors
## Filtering anchors
##  Retained 2412 anchors
# Create the integrated data assay
scRNA_combined <- IntegrateData(anchorset=brca_anchors)
## Merging dataset 4 into 3
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
## Merging dataset 1 into 2
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
## Merging dataset 2 1 into 3 4
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data

Quality control

# Use the integrated assay as default
DefaultAssay(scRNA_combined) <- "RNA"

# Creating a new slot with the percentage of counts on mitochondrial reads


scRNA_combined[["percent.mt"]] <- PercentageFeatureSet(scRNA_combined, pattern="MT-")

# Violin plots of specific features
VlnPlot(scRNA_combined, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3,
        split.by="orig.ident")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

We will remove cells with a percentage of mitochondrial reads too high (>5%) and with an abnormally high or low number of feature counts, that can suggest empty droplets or multiplets in the microfluidic system, respectively.

scRNA_combined <- subset(scRNA_combined, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 10)

VlnPlot(scRNA_combined, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3,
        split.by="orig.ident")

Visualization and clustering

Check integration

# Change the default to the integrated dataset
DefaultAssay(scRNA_combined) <- "integrated"

scRNA_combined <- ScaleData(scRNA_combined, verbose = FALSE)
scRNA_combined <- RunPCA(scRNA_combined, npcs = 30, verbose = FALSE)
scRNA_combined <- RunUMAP(scRNA_combined, reduction = "pca", dims = 1:30)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 22:19:16 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:19:16 Read 3006 rows and found 30 numeric columns
## 22:19:16 Using Annoy for neighbor search, n_neighbors = 30
## 22:19:16 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:19:17 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c705eea473b
## 22:19:17 Searching Annoy index using 1 thread, search_k = 3000
## 22:19:18 Annoy recall = 100%
## 22:19:19 Commencing smooth kNN distance calibration using 1 thread
## 22:19:21 Initializing from normalized Laplacian + noise
## 22:19:22 Commencing optimization for 500 epochs, with 129862 positive edges
## 22:19:34 Optimization finished
scRNA_combined <- FindNeighbors(scRNA_combined, reduction = "pca", dims = 1:30)
## Computing nearest neighbor graph
## Computing SNN
scRNA_combined <- FindClusters(scRNA_combined, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 3006
## Number of edges: 151626
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7945
## Number of communities: 5
## Elapsed time: 0 seconds
# Visualization
p1 <- DimPlot(scRNA_combined, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(scRNA_combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

DimPlot(scRNA_combined, reduction = "umap", split.by = "orig.ident")

After normalization and integration, the UMAPs of all subsets look quite similar, so we can assume integration is correct and we can proceed forward. Once this has been verified, we can go back a few steps and refine our analysis

Determine the dimensionality of the dataset

scRNA_combined <- JackStraw(scRNA_combined, num.replicate = 100)
scRNA_combined <- ScoreJackStraw(scRNA_combined, dims=1:20)

JackStrawPlot(scRNA_combined, dims=1:20)
## Warning: Removed 28000 rows containing missing values (geom_point).

ElbowPlot(scRNA_combined)

Based on these two graphs, we may consider using 13-15 PCs.

scRNA_combined <- FindNeighbors(scRNA_combined, reduction = "pca", dims = 1:13)
## Computing nearest neighbor graph
## Computing SNN
scRNA_combined <- FindClusters(scRNA_combined, resolution = 0.6) # Try different levels of resolution and chose the one that results in cluster with most biological sense1
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 3006
## Number of edges: 112761
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7863
## Number of communities: 6
## Elapsed time: 0 seconds
scRNA_combined <- RunUMAP(scRNA_combined, dims = 1:13)
## 22:22:57 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:22:57 Read 3006 rows and found 13 numeric columns
## 22:22:57 Using Annoy for neighbor search, n_neighbors = 30
## 22:22:57 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:22:58 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c7055085643
## 22:22:58 Searching Annoy index using 1 thread, search_k = 3000
## 22:23:00 Annoy recall = 100%
## 22:23:01 Commencing smooth kNN distance calibration using 1 thread
## 22:23:03 Initializing from normalized Laplacian + noise
## 22:23:03 Commencing optimization for 500 epochs, with 123068 positive edges
## 22:23:17 Optimization finished
DimPlot(scRNA_combined, reduction="umap")

DimPlot(scRNA_combined, reduction = "umap", split.by = "orig.ident")

Finding marker genes

hg19_markers <- FindAllMarkers(scRNA_combined, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
hg19_markers %>%
    group_by(cluster) %>%
    slice_max(n = 2, order_by = avg_log2FC)
## # A tibble: 12 x 7
## # Groups:   cluster [6]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene    
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
##  1 3.69e- 69      0.815 0.926 0.82  7.38e- 66 0       CLEC3A  
##  2 3.34e- 72      0.767 0.939 0.843 6.69e- 69 0       AGR2    
##  3 1.15e- 79      0.879 0.997 0.931 2.30e- 76 1       NDRG1   
##  4 3.76e- 30      0.500 0.985 0.889 7.51e- 27 1       IGFBP5  
##  5 2.28e-  6      1.55  0.716 0.73  4.57e-  3 2       SCGB2A2 
##  6 3.31e- 88      1.11  0.973 0.896 6.62e- 85 2       IGFBP5  
##  7 1.53e-161      2.03  0.969 0.7   3.06e-158 3       HIST1H4C
##  8 2.69e-218      1.93  0.971 0.632 5.38e-215 3       KIAA0101
##  9 5.22e-142      3.11  0.951 0.561 1.04e-138 4       KRT16   
## 10 5.35e-135      2.87  0.904 0.447 1.07e-131 4       KRT6B   
## 11 6.67e- 90      2.70  0.98  0.71  1.33e- 86 5       UBE2C   
## 12 8.13e-113      2.28  0.995 0.686 1.63e-109 5       PTTG1
hg19_markers %>%
    group_by(cluster) %>%
    top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(scRNA_combined, features = top10$gene) + NoLegend()

Functional analysis

Gene set enrichment analysis

MSigDB_df <- msigdbr(species = "Homo sapiens")
fgsea_sets<- MSigDB_df %>% split(x = .$gene_symbol, f = .$gs_name)



plot_GSEA <- function(cluster_number){
  genes <- hg19_markers %>% filter(cluster==as.character(cluster_number)) %>%
  arrange(desc(p_val_adj)) %>% 
  dplyr::select(gene, avg_log2FC)
  
  # Create dataframe with genes
  ranks <- deframe(genes)
  
  # Perform GSEA analysis
  fgseaRes <- fgsea(fgsea_sets,
                    stats = ranks,
                    minSize=10,
                    maxSize=500,
                    nperm=1000000)
  
  # Plot
  ggplot(fgseaRes %>% filter(padj < 0.01) %>% head(n= 20), aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill= NES < 7.5)) +
  coord_flip() +
  labs(x="Pathway",
       y="Normalized Enrichment Score",
       title=paste0("Hallmark pathways NES from GSEA for cluster",as.character(cluster_number)))
}

for (cluster in 0:5){plot_GSEA(cluster)}

Gene ontology

plot_GO <- function(cluster_number){
  genes <- hg19_markers %>% filter(cluster==as.character(cluster_number))
  
  # Calculate enrichment
  enrichment <- enrichGO(gene = genes$gene,
             OrgDb  = org.Hs.eg.db,
             keyType = "SYMBOL",
             ont  = "BP",
             pAdjustMethod = "BH",
             pvalueCutoff  = 0.01,
             qvalueCutoff  = 0.05)
  enrichment <- clusterProfiler::simplify(enrichment, cutoff=0.7, by="p.adjust", select_fun=min)
  
  GO_ggdata <- enrichment %>%
   as_data_frame() %>%
   arrange(Count)
  GO_ggdata$Description <- factor(GO_ggdata$Description, levels = GO_ggdata$Description)

ggplot(GO_ggdata, aes(x = Description, y = Count, fill = p.adjust)) +
 geom_bar(stat = "identity") +
 scale_colour_viridis_d(begin=0,end=1) +
 coord_flip() +
 ylab("Number of genes") +
 xlab("GO Terms") +
 theme(axis.text.y = element_text(size=10))
}

for (cluster in 0:5){print(plot_GO(cluster))}
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.

Automated cluster annotation

Method 1- scCATCH

hg19_geneinfo <- rev_gene(data = scRNA_combined[['integrated']]@data,
                          data_type = "data",
                          species="Human",
                          geneinfo = geneinfo)

hg19_ann_1 <- createscCATCH(data = scRNA_combined[['integrated']]@data, 
                          cluster = as.character(Idents(scRNA_combined)))

hg19_ann_1 <- findmarkergene(object = hg19_ann_1, species = "Human", marker = cellmatch, tissue = "Breast", cancer="Breast Cancer", use_method = "2")
## There are 535 potential marker genes in CellMatch database for Human Breast Cancer on Breast.
hg19_ann_1 <- findcelltype(hg19_ann_1)
hg19_ann_1@celltype
##   cluster                                                  cluster_marker
## 1       1                              FN1, IGFBP5, TACSTD2, PMEPA1, CD44
## 2       2                                                          IGFBP5
## 3       3                                      CENPF, BIRC5, HELLS, BLVRB
## 4       0                                                            BST2
## 5       5                                                    CENPF, BIRC5
## 6       4 G0S2, FN1, ITGA6, TACSTD2, LIMA1, PMEPA1, DNAJB1, ZFP36L2, CD44
##          cell_type celltype_score celltype_related_marker
## 1 Cancer Stem Cell           0.98                    CD44
## 2      Killer Cell           0.50                  IGFBP5
## 3    Helper T Cell           0.61     CENPF, BIRC5, HELLS
## 4           T Cell           0.50                    BST2
## 5    Helper T Cell           0.58            CENPF, BIRC5
## 6 Cancer Stem Cell           0.98             CD44, ITGA6
##                                                                                                                                                                                                                                                                                                                                                                                                   PMID
## 1 24099815, 28986882, 27840965, 23799994, 23053657, 23768049, 28719381, 28579829, 28554844, 28399903, 28166194, 28148288, 28036386, 28012234, 27768764, 27630305, 27458252, 27133070, 27115469, 26893363, 26170011, 25990436, 24596390, 24297508, 24171482, 24144739, 24055090, 23543868, 23112116, 23056395, 23036368, 22761812, 22464177, 21979823, 21802218, 21092249, 22459349, 27640754, 25940879
## 2                                                                                                                                                                                                                                                                                                                                                                                             28474673
## 3                                                                                                                                                                                                                                                                                                                                                                                             28474673
## 4                                                                                                                                                                                                                                                                                                                                                                                             28474673
## 5                                                                                                                                                                                                                                                                                                                                                                                             28474673
## 6 24099815, 28986882, 27840965, 23799994, 23053657, 23768049, 28719381, 28579829, 28554844, 28399903, 28166194, 28148288, 28036386, 28012234, 27768764, 27630305, 27458252, 27133070, 27115469, 26893363, 26170011, 25990436, 24596390, 24297508, 24171482, 24144739, 24055090, 23543868, 23112116, 23056395, 23036368, 22761812, 22464177, 21979823, 21802218, 21092249, 22459349, 27640754, 25940879

Method 2- SingleR

hpca.se <- HumanPrimaryCellAtlasData()
## Warning: 'HumanPrimaryCellAtlasData' is deprecated.
## Use 'celldex::HumanPrimaryCellAtlasData' instead.
## See help("Deprecated")
## snapshotDate(): 2020-10-27
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
hpca.se
## class: SummarizedExperiment 
## dim: 19363 713 
## metadata(0):
## assays(1): logcounts
## rownames(19363): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
## rowData names(0):
## colnames(713): GSM112490 GSM112491 ... GSM92233 GSM92234
## colData names(3): label.main label.fine label.ont
hg19_ann_2 <- SingleR(test = scRNA_combined[['integrated']]@data, 
                      ref = hpca.se, 
                      assay.type.test=1,
                      labels = hpca.se$label.main)

hg19_ann_2
## DataFrame with 3006 rows and 5 columns
##                                              scores      first.labels
##                                            <matrix>       <character>
## AAACCTGGTTCCACTC-1_1 0.284388:0.263078:0.208514:...  Epithelial_cells
## AACCGCGCACCTTGTC-1_1 0.219088:0.219750:0.185685:...  Epithelial_cells
## AACCGCGTCCAAGCCG-1_1 0.288097:0.332301:0.294157:... Endothelial_cells
## AACGTTGGTTCAGTAC-1_1 0.312886:0.359925:0.326496:...        BM & Prog.
## AACGTTGTCCTGCTTG-1_1 0.254901:0.215151:0.166806:...  Epithelial_cells
## ...                                             ...               ...
## TTCTTAGCATTCTTAC-1_4 0.128125:0.109205:0.078661:...       Hepatocytes
## TTGACTTGTATGCTTG-1_4 0.171581:0.194234:0.169912:...       Hepatocytes
## TTGACTTTCCCACTTG-1_4 0.175745:0.202312:0.173485:...  Epithelial_cells
## TTGGCAAGTCTTGTCC-1_4 0.192326:0.156040:0.144870:...  Epithelial_cells
## TTTACTGGTAGGCATG-1_4 0.140242:0.201939:0.168875:...  Epithelial_cells
##                            tuning.scores            labels     pruned.labels
##                              <DataFrame>       <character>       <character>
## AAACCTGGTTCCACTC-1_1 0.223729: 0.2155909  Epithelial_cells  Epithelial_cells
## AACCGCGCACCTTGTC-1_1 0.244705:-0.0810895  Epithelial_cells  Epithelial_cells
## AACCGCGTCCAAGCCG-1_1 0.183081: 0.1802609 Endothelial_cells Endothelial_cells
## AACGTTGGTTCAGTAC-1_1 0.297462: 0.2850939        BM & Prog.        BM & Prog.
## AACGTTGTCCTGCTTG-1_1 0.220434: 0.1080830  Epithelial_cells  Epithelial_cells
## ...                                  ...               ...               ...
## TTCTTAGCATTCTTAC-1_4 0.1070612:0.0436891  Epithelial_cells  Epithelial_cells
## TTGACTTGTATGCTTG-1_4 0.0747311:0.0635165       Hepatocytes       Hepatocytes
## TTGACTTTCCCACTTG-1_4 0.3865267:0.3086367  Epithelial_cells  Epithelial_cells
## TTGGCAAGTCTTGTCC-1_4 0.2565691:0.1856667  Epithelial_cells  Epithelial_cells
## TTTACTGGTAGGCATG-1_4 0.1900828:0.1195525  Epithelial_cells  Epithelial_cells

Analysis of mouse cells

Data integration

# Normalize and identify variable features
scRNA_files_mm10 <- lapply(X = scRNA_files_mm10, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# Select features repeatedly variable between the 4 datasets
features_mm10 <- SelectIntegrationFeatures(object.list = scRNA_files_mm10)

# Select anchors for integration
brca_anchors_mm10 <- FindIntegrationAnchors(object.list = scRNA_files_mm10, anchor.features = features_mm10)
## Warning in CheckDuplicateCellNames(object.list = object.list): Some cell names
## are duplicated across objects provided. Renaming to enforce unique cell names.
## Scaling features for provided objects
## Finding all pairwise anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2294 anchors
## Filtering anchors
##  Retained 2030 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2088 anchors
## Filtering anchors
##  Retained 1826 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2833 anchors
## Filtering anchors
##  Retained 2294 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2231 anchors
## Filtering anchors
##  Retained 2108 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2983 anchors
## Filtering anchors
##  Retained 2716 anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 2769 anchors
## Filtering anchors
##  Retained 2515 anchors
# Create the integrated data assay
scRNA_combined_mm10 <- IntegrateData(anchorset=brca_anchors_mm10)
## Merging dataset 1 into 4
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
## Merging dataset 2 into 4 1
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
## Merging dataset 3 into 4 1 2
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data

Quality control

# Use the integrated assay as default
DefaultAssay(scRNA_combined_mm10) <- "RNA"

# Creating a new slot with the percentage of counts on mitochondrial reads


scRNA_combined_mm10[["percent.mt"]] <- PercentageFeatureSet(scRNA_combined_mm10, pattern="mt-")

# Violin plots of specific features
VlnPlot(scRNA_combined_mm10, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3,
        split.by="orig.ident")

We will remove cells with a percentage of mitochondrial reads too high (>5%) and with an abnormally high or low number of feature counts, that can suggest empty droplets or multiplets in the microfluidic system, respectively.

scRNA_combined_mm10 <- subset(scRNA_combined_mm10, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 10)

VlnPlot(scRNA_combined_mm10, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3,
        split.by="orig.ident")

Visualization and clustering

Check integration

# Change the default to the integrated dataset
DefaultAssay(scRNA_combined_mm10) <- "integrated"

scRNA_combined_mm10 <- ScaleData(scRNA_combined_mm10, verbose = FALSE)
scRNA_combined_mm10 <- RunPCA(scRNA_combined_mm10, npcs = 30, verbose = FALSE)
scRNA_combined_mm10 <- RunUMAP(scRNA_combined_mm10, reduction = "pca", dims = 1:30)
## 22:31:28 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:31:28 Read 4017 rows and found 30 numeric columns
## 22:31:28 Using Annoy for neighbor search, n_neighbors = 30
## 22:31:28 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:31:29 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c70633baa4
## 22:31:29 Searching Annoy index using 1 thread, search_k = 3000
## 22:31:31 Annoy recall = 100%
## 22:31:33 Commencing smooth kNN distance calibration using 1 thread
## 22:31:35 Initializing from normalized Laplacian + noise
## 22:31:36 Commencing optimization for 500 epochs, with 167364 positive edges
## 22:31:56 Optimization finished
scRNA_combined_mm10 <- FindNeighbors(scRNA_combined_mm10, reduction = "pca", dims = 1:30)
## Computing nearest neighbor graph
## Computing SNN
scRNA_combined_mm10 <- FindClusters(scRNA_combined_mm10, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4017
## Number of edges: 155514
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9021
## Number of communities: 14
## Elapsed time: 0 seconds
# Visualization
p1 <- DimPlot(scRNA_combined_mm10, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(scRNA_combined_mm10, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

DimPlot(scRNA_combined_mm10, reduction = "umap", split.by = "orig.ident")

After normalization and integration, the UMAPs of all subsets look quite similar, so we can assume integration is correct and we can proceed forward. Once this has been verified, we can go back a few steps and refine our analysis

Determine the dimensionality of the dataset

scRNA_combined_mm10 <- JackStraw(scRNA_combined_mm10, num.replicate = 100)
scRNA_combined_mm10 <- ScoreJackStraw(scRNA_combined_mm10, dims=1:20)

JackStrawPlot(scRNA_combined_mm10, dims=1:20)
## Warning: Removed 28000 rows containing missing values (geom_point).

ElbowPlot(scRNA_combined_mm10)

Based on these two graphs, we may consider using 10 PCs.

scRNA_combined_mm10 <- FindNeighbors(scRNA_combined_mm10, reduction = "pca", dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
scRNA_combined_mm10 <- FindClusters(scRNA_combined_mm10, resolution = 0.6) # Try different levels of resolution and chose the one that results in cluster with most biological sense1
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4017
## Number of edges: 128730
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8992
## Number of communities: 12
## Elapsed time: 0 seconds
scRNA_combined_mm10 <- RunUMAP(scRNA_combined_mm10, dims = 1:10)
## 22:36:02 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:36:02 Read 4017 rows and found 10 numeric columns
## 22:36:02 Using Annoy for neighbor search, n_neighbors = 30
## 22:36:02 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:36:03 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c7063c14705
## 22:36:03 Searching Annoy index using 1 thread, search_k = 3000
## 22:36:05 Annoy recall = 100%
## 22:36:07 Commencing smooth kNN distance calibration using 1 thread
## 22:36:10 Initializing from normalized Laplacian + noise
## 22:36:10 Commencing optimization for 500 epochs, with 160052 positive edges
## 22:36:29 Optimization finished
DimPlot(scRNA_combined_mm10, reduction="umap")

DimPlot(scRNA_combined_mm10, reduction = "umap", split.by = "orig.ident")

Finding marker genes

mm10_markers <- FindAllMarkers(scRNA_combined_mm10, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 11
mm10_markers %>%
    group_by(cluster) %>%
    slice_max(n = 2, order_by = avg_log2FC)
## # A tibble: 24 x 7
## # Groups:   cluster [12]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr> 
##  1 4.20e-194       2.98 0.868 0.692 8.39e-191 0       Gzme  
##  2 3.12e-183       2.86 0.96  0.748 6.24e-180 0       Spp1  
##  3 1.17e-304       3.83 0.994 0.742 2.35e-301 1       C1qa  
##  4 0               3.61 1     0.725 0         1       Ctss  
##  5 6.73e-127       2.38 0.871 0.597 1.35e-123 2       Mmp3  
##  6 2.57e-223       2.10 0.945 0.547 5.14e-220 2       Mfap4 
##  7 6.65e-190       2.73 0.976 0.647 1.33e-186 3       Apod  
##  8 3.04e-184       2.55 0.988 0.797 6.08e-181 3       Gsn   
##  9 6.90e-169       3.22 0.997 0.696 1.38e-165 4       Clec3b
## 10 1.34e-169       3.21 0.983 0.409 2.67e-166 4       Sbsn  
## # ... with 14 more rows
mm10_markers %>%
    group_by(cluster) %>%
    top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(scRNA_combined_mm10, features = top10$gene) + NoLegend()

Functional analysis

Gene set enrichment analysis

MSigDB_df <- msigdbr(species = "Homo sapiens")
fgsea_sets<- MSigDB_df %>% split(x = .$gene_symbol, f = .$gs_name)



plot_GSEA <- function(cluster_number){
  genes <- mm10_markers %>% filter(cluster==as.character(cluster_number)) %>%
  arrange(desc(p_val_adj)) %>% 
  dplyr::select(gene, avg_log2FC)
  
  # Create dataframe with genes
  ranks <- deframe(genes)
  
  # Perform GSEA analysis
  fgseaRes <- fgsea(fgsea_sets,
                    stats = ranks,
                    minSize=10,
                    maxSize=500,
                    nperm=1000000)
  
  # Plot
  ggplot(fgseaRes %>% filter(padj < 0.01) %>% head(n= 20), aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill= NES < 7.5)) +
  coord_flip() +
  labs(x="Pathway",
       y="Normalized Enrichment Score",
       title=paste0("Hallmark pathways NES from GSEA for cluster",as.character(cluster_number)))
}

plot_GSEA(0)

for (cluster in 0:5){plot_GSEA(cluster)}

Gene ontology

plot_GO_mm10 <- function(cluster_number){
  genes <- mm10_markers %>% filter(cluster==as.character(cluster_number))
  
  # Calculate enrichment
  enrichment <- enrichGO(gene = genes$gene,
             OrgDb  = org.Mm.eg.db,
             keyType = "SYMBOL",
             ont  = "BP",
             pAdjustMethod = "BH",
             pvalueCutoff  = 0.01,
             qvalueCutoff  = 0.05)
  enrichment <- clusterProfiler::simplify(enrichment, cutoff=0.7, by="p.adjust", select_fun=min)
  
  GO_ggdata <- enrichment %>%
   as_data_frame() %>%
   arrange(Count)
  GO_ggdata$Description <- factor(GO_ggdata$Description, levels = GO_ggdata$Description)

ggplot(GO_ggdata, aes(x = Description, y = Count, fill = p.adjust)) +
 geom_bar(stat = "identity") +
 scale_colour_viridis_d(begin=0,end=1) +
 coord_flip() +
 ylab("Number of genes") +
 xlab("GO Terms") +
 theme(axis.text.y = element_text(size=10))
}

for (cluster in 0:5){print(plot_GO_mm10(cluster))}

Automated cluster annotation

Method 1- scCATCH

MouseRNAseqData(ensembl = FALSE, cell.ont = c("all", "nonna", "none"))
## snapshotDate(): 2020-10-27
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
## class: SummarizedExperiment 
## dim: 21214 358 
## metadata(0):
## assays(1): logcounts
## rownames(21214): Xkr4 Rp1 ... LOC100039574 LOC100039753
## rowData names(0):
## colnames(358): ERR525589Aligned ERR525592Aligned ... SRR1044043Aligned
##   SRR1044044Aligned
## colData names(3): label.main label.fine label.ont
mm10_geneinfo <- rev_gene(data = scRNA_combined_mm10[['integrated']]@data,
                          data_type = "data",
                          species="Mouse",
                          geneinfo = geneinfo)

mm10_ann_1 <- createscCATCH(data = scRNA_combined_mm10[['integrated']]@data, 
                          cluster = as.character(Idents(scRNA_combined_mm10)))

mm10_ann_1 <- findmarkergene(object = mm10_ann_1, species = "Mouse", marker = cellmatch, tissue="Mammary gland", use_method = "2")
## There are 380 potential marker genes in CellMatch database for Mouse on Mammary gland.
mm10_ann_1 <- findcelltype(mm10_ann_1)
mm10_ann_1@celltype
##    cluster
## 1        4
## 2        5
## 3        0
## 4        7
## 5        1
## 6        3
## 7        6
## 8        2
## 9        8
## 10       9
## 11      11
## 12      10
##                                                                                                                                                                                                                                                                                                                                                                                             cluster_marker
## 1                                                                              Igfbp5, Apod, Mt1, Clec3b, Gsn, Aldh1a3, Dcn, Plet1, Serpinb1a, Igfbp4, Efemp1, Procr, Igfbp6, Rarres2, Fn1, Htra3, S100a4, Ogn, Dpt, Ly6c1, Col3a1, Lum, Gadd45g, Ctsl, Mfap5, Serping1, Glul, Gpx3, Ly6a, Dpep1, Serpinf1, Ctsk, Pcolce, Loxl1, Fstl1, Phlda1, Plpp3, S100a6, Aebp1, Pdpn, C1s1, Mmp2, Lsp1, Anxa1, Anxa2
## 2                                                                                                                                                                                        Fabp4, Ptn, Gng11, Lrg1, Col4a2, Atf3, Procr, Col4a1, Ly6c1, Id3, Rbp1, Fscn1, Crip2, Sparc, Csrp1, Tspan13, Ly6a, Cdkn1a, Tsc22d1, Kit, Itga6, Cd200, Gata2, Cnn3, Slc12a2, Igfbp7, S100a6, Serpinh1, Vim, Anxa2
## 3                                               Actg2, Gng11, Acta2, Igfbp2, Tagln, Ctsd, Myl9, Cnn1, Col4a2, Postn, Tpm2, Igfbp6, Tpm1, Col6a3, Col1a1, Col1a2, Fn1, Col3a1, Rbp1, Fscn1, Crip2, Col6a1, Mfap5, Sparc, Csrp1, Mfge8, Col6a2, Gpx3, Bgn, Lgals1, Tsc22d1, Serpinf1, Col5a2, Cd200, Cnn3, Pcolce, Loxl1, Fstl1, Mmp14, Igfbp7, S100a6, Rcn3, Serpinh1, Aebp1, Vim, Mmp2, Lsp1, Anxa1, Anxa2
## 4                                                                                                                                                                                                                                   Ccl5, Hp, H2-Ab1, Lyz2, Ccl6, H2-Aa, Ctss, Cd14, Atf3, Cd68, H2-Eb1, AW112010, S100a4, Cd24a, Fcer1g, Cfp, Trf, Aif1, Tyrobp, Ms4a6c, Tspan13, Cdkn1a, Ctsc, Tgm2, Grn
## 5                                                                                                                                                                                                                                               Ccl5, H2-Ab1, Lyz2, H2-Aa, Ctss, C1qc, Ctsd, Cd14, Cd68, H2-Eb1, AW112010, Fcer1g, Trf, Aif1, Tyrobp, Ms4a6c, Tspan13, Cdkn1a, Ctsc, Tgm2, Ctsb, Lgmn, Grn
## 6    Cxcl14, Apod, Mt1, Lcn2, Lpl, Clec3b, Gsn, Dcn, Serpinb1a, Igfbp4, Efemp1, Procr, Igfbp6, Col4a1, Rarres2, Col1a1, Col1a2, Fn1, Htra3, Ogn, Dpt, Ly6c1, Col3a1, Id3, Lum, Gadd45g, Ctsl, Col6a1, Mfap5, Sparc, Serping1, Glul, Col6a2, Gpx3, Ly6a, Dpep1, Bgn, Serpinf1, Col5a2, Ctsk, Pcolce, Mfap2, Loxl1, Fstl1, Phlda1, Plpp3, S100a6, Rcn3, Serpinh1, Aebp1, Pdpn, Dbi, C1s1, Mmp2, Anxa1, Anxa2
## 7                                                                                                                                  Actg2, Gng11, Acta2, Tagln, Mylk, Myl9, Cnn1, Col4a2, Postn, Tpm2, Tpm1, Col4a1, Rarres2, Col6a3, Col1a1, Col1a2, Fn1, S100a4, Col3a1, Crip2, Sparc, Csrp1, Mfge8, Bgn, Lgals1, Col5a2, Pcolce, Fstl1, Igfbp7, Cystm1, S100a6, Rcn3, Serpinh1, Aebp1, Vim, Anxa1, Anxa2
## 8  Cxcl14, Apod, Mt1, Dcn, Postn, Igfbp4, Igfbp6, Tpm1, Itm2a, Rarres2, Col6a3, Col1a1, Col1a2, Fn1, Ogn, Dpt, Col3a1, Rbp1, Lum, Gadd45g, Ctsl, Col6a1, Mfap5, Sparc, Serping1, Glul, Mfge8, Col6a2, Gpx3, Dpep1, Bgn, Lgals1, Tsc22d1, Serpinf1, Col5a2, Ctsk, Cnn3, Pcolce, Mfap2, Loxl1, Fstl1, Mmp14, Phlda1, S100a6, Rcn3, Serpinh1, Aebp1, Pdpn, Dbi, Vim, C1s1, Tmem176a, Mmp2, Lsp1, Anxa1, Anxa2
## 9                                                                                                                                                                                                                                                                                                                                                                            Ccl5, Nkg7, AW112010, Tspan13
## 10                                                                                                                                                                                                          H2-Ab1, Lyz2, Ccl6, H2-Aa, Pf4, Mt1, Ctss, C1qc, Ctsd, Cd14, Atf3, Cd68, H2-Eb1, AW112010, Fcer1g, Cfp, Trf, Aif1, Tyrobp, Ms4a6c, Fxyd2, Cdkn1a, Ctsc, Ctsb, Cystm1, Dbi, Lgmn, Tmem176a, Grn
## 11                                                                                                                                                                                                                                                H2-Ab1, Lyz2, H2-Aa, Ctss, C1qc, Ctsd, Cd14, Cd68, H2-Eb1, AW112010, Fcer1g, Cfp, Trf, Tuba1b, Aif1, Tyrobp, Ms4a6c, Cdkn1a, Ctsc, Tgm2, Ctsb, Lgmn, Grn
## 12                                                                                                            Actg2, Gng11, Acta2, Tagln, Myl9, Dcn, Postn, Tpm2, Tpm1, Col4a1, Col6a3, Col1a1, Col1a2, Fn1, S100a4, Col3a1, Tuba1b, Lum, Crip2, Col6a1, Sparc, Csrp1, Mfge8, Col6a2, Bgn, Lgals1, Serpinf1, Col5a2, Cnn3, Pcolce, Loxl1, Fstl1, Igfbp7, Phlda1, S100a6, Rcn3, Serpinh1, Vim, Anxa1, Anxa2
##                    cell_type celltype_score
## 1    Luminal Progenitor Cell           0.83
## 2    Luminal Progenitor Cell           0.88
## 3         Myoepithelial Cell           0.79
## 4    Luminal Progenitor Cell           0.79
## 5    Luminal Progenitor Cell           0.73
## 6    Luminal Progenitor Cell           0.81
## 7         Myoepithelial Cell           0.79
## 8    Luminal Progenitor Cell           0.76
## 9  Killer Cell, Luminal Cell     0.58, 0.58
## 10   Luminal Progenitor Cell           0.73
## 11   Luminal Progenitor Cell           0.73
## 12        Myoepithelial Cell           0.78
##                                                                      celltype_related_marker
## 1                                                Aldh1a3, Plet1, Igfbp5, Glul, Anxa2, Phlda1
## 2                                                    Kit, Lrg1, Slc12a2, Ptn, Tsc22d1, Anxa2
## 3  Acta2, Tagln, Actg2, Igfbp2, Cnn1, Myl9, Tpm2, Tpm1, Csrp1, Rbp1, Mfge8, Col4a2, Serpinh1
## 4                                                                                 Cd14, Tgm2
## 5                                                                                 Cd14, Tgm2
## 6                                                             Glul, Lcn2, Dbi, Anxa2, Phlda1
## 7  Acta2, Tagln, Actg2, Cnn1, Mylk, Myl9, Tpm2, Tpm1, Csrp1, Mfge8, Col4a1, Col4a2, Serpinh1
## 8                                                   Mfge8, Glul, Tsc22d1, Dbi, Anxa2, Phlda1
## 9                                                              AW112010, Ccl5, Tspan13, Nkg7
## 10                                                                                 Cd14, Dbi
## 11                                                                                Cd14, Tgm2
## 12                     Acta2, Tagln, Actg2, Myl9, Tpm2, Tpm1, Csrp1, Mfge8, Col4a1, Serpinh1
##                            PMID
## 1            29225342, 29775597
## 2  21132007, 29225342, 29775597
## 3            29225342, 29775597
## 4            29225342, 29775597
## 5            29225342, 29775597
## 6                      29775597
## 7            29225342, 29775597
## 8                      29775597
## 9                      29775597
## 10           29225342, 29775597
## 11           29225342, 29775597
## 12           29225342, 29775597

Method 2- SingleR

mRNAseq.se <- MouseRNAseqData()
## snapshotDate(): 2020-10-27
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
mRNAseq.se
## class: SummarizedExperiment 
## dim: 21214 358 
## metadata(0):
## assays(1): logcounts
## rownames(21214): Xkr4 Rp1 ... LOC100039574 LOC100039753
## rowData names(0):
## colnames(358): ERR525589Aligned ERR525592Aligned ... SRR1044043Aligned
##   SRR1044044Aligned
## colData names(3): label.main label.fine label.ont
mm10_ann_2 <- SingleR(test = scRNA_combined_mm10[['integrated']]@data, 
                      ref = mRNAseq.se, 
                      assay.type.test=1,
                      labels = mRNAseq.se$label.main)

mm10_ann_2
## DataFrame with 4017 rows and 5 columns
##                                                    scores      first.labels
##                                                  <matrix>       <character>
## AAACCTGCATCACCCT-1_1 0.2749191: 0.0953161:-0.00695552:...       Fibroblasts
## AAACCTGCATTGGGCC-1_1 0.1399648: 0.0770076: 0.09786590:... Endothelial cells
## AAACCTGGTGGTTTCA-1_1 0.3655346: 0.2011414: 0.07613095:...       Fibroblasts
## AAAGCAAAGGATCGCA-1_1 0.0561536:-0.0161654: 0.16629786:...         Monocytes
## AAAGCAAAGTCATGCT-1_1 0.0779979: 0.0166078: 0.10401372:...         Monocytes
## ...                                                   ...               ...
## TTTGGTTGTGACGGTA-1_4   0.37394650:0.1605742:0.0708835:...       Fibroblasts
## TTTGGTTGTTGTCTTT-1_4   0.25604361:0.1082203:0.1364289:...       Fibroblasts
## TTTGGTTTCTTTAGGG-1_4   0.30681930:0.1277484:0.0688820:...       Fibroblasts
## TTTGTCAAGAGTGAGA-1_4   0.00969858:0.0060175:0.3195598:...           T cells
## TTTGTCAAGCTTATCG-1_4   0.35637833:0.1532604:0.0670055:...       Fibroblasts
##                            tuning.scores            labels     pruned.labels
##                              <DataFrame>       <character>       <character>
## AAACCTGCATCACCCT-1_1   0.326946:0.274919       Fibroblasts       Fibroblasts
## AAACCTGCATTGGGCC-1_1   0.218958:0.144200 Endothelial cells Endothelial cells
## AAACCTGGTGGTTTCA-1_1   0.432832:0.365535       Fibroblasts       Fibroblasts
## AAAGCAAAGGATCGCA-1_1   0.366443:0.228990         Monocytes         Monocytes
## AAAGCAAAGTCATGCT-1_1   0.284693:0.222767         Monocytes         Monocytes
## ...                                  ...               ...               ...
## TTTGGTTGTGACGGTA-1_4 0.385029:0.07043289       Fibroblasts       Fibroblasts
## TTTGGTTGTTGTCTTT-1_4 0.409715:0.00152361       Fibroblasts       Fibroblasts
## TTTGGTTTCTTTAGGG-1_4 0.361933:0.30681930       Fibroblasts       Fibroblasts
## TTTGTCAAGAGTGAGA-1_4 0.353053:0.28400529           T cells           T cells
## TTTGTCAAGCTTATCG-1_4 0.437798:0.35637833       Fibroblasts       Fibroblasts

Session info

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.1252    
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] celldex_1.0.0               org.Mm.eg.db_3.12.0        
##  [3] org.Hs.eg.db_3.12.0         AnnotationDbi_1.52.0       
##  [5] clusterProfiler_3.18.1      fgsea_1.16.0               
##  [7] tibble_3.1.6                msigdbr_7.5.1              
##  [9] SingleR_1.4.1               SummarizedExperiment_1.20.0
## [11] Biobase_2.50.0              GenomicRanges_1.42.0       
## [13] GenomeInfoDb_1.26.7         IRanges_2.24.1             
## [15] S4Vectors_0.28.1            BiocGenerics_0.36.1        
## [17] MatrixGenerics_1.2.1        matrixStats_0.62.0         
## [19] BiocManager_1.30.17         ggplot2_3.3.5              
## [21] scCATCH_3.2                 purrr_0.3.4                
## [23] sp_1.4-7                    SeuratObject_4.1.0         
## [25] Seurat_4.1.1                stringr_1.4.0              
## [27] dplyr_1.0.8                 R.utils_2.11.0             
## [29] R.oo_1.24.0                 R.methodsS3_1.8.1          
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2                    reticulate_1.25              
##   [3] tidyselect_1.1.2              RSQLite_2.2.12               
##   [5] htmlwidgets_1.5.4             grid_4.0.4                   
##   [7] BiocParallel_1.24.1           Rtsne_0.16                   
##   [9] scatterpie_0.1.7              munsell_0.5.0                
##  [11] codetools_0.2-18              ica_1.0-2                    
##  [13] future_1.25.0                 miniUI_0.1.1.1               
##  [15] withr_2.5.0                   spatstat.random_2.2-0        
##  [17] colorspace_2.0-3              GOSemSim_2.16.1              
##  [19] progressr_0.10.0              highr_0.9                    
##  [21] knitr_1.31                    rstudioapi_0.13              
##  [23] ROCR_1.0-11                   tensor_1.5                   
##  [25] DOSE_3.16.0                   listenv_0.8.0                
##  [27] labeling_0.4.2                GenomeInfoDbData_1.2.4       
##  [29] polyclip_1.10-0               farver_2.1.0                 
##  [31] bit64_4.0.5                   downloader_0.4               
##  [33] parallelly_1.31.1             vctrs_0.4.1                  
##  [35] generics_0.1.2                xfun_0.21                    
##  [37] BiocFileCache_1.14.0          R6_2.5.1                     
##  [39] graphlayouts_0.8.0            rsvd_1.0.5                   
##  [41] bitops_1.0-7                  spatstat.utils_2.3-1         
##  [43] cachem_1.0.4                  DelayedArray_0.16.3          
##  [45] assertthat_0.2.1              promises_1.2.0.1             
##  [47] scales_1.2.0                  ggraph_2.0.5                 
##  [49] enrichplot_1.10.2             rgeos_0.5-9                  
##  [51] gtable_0.3.0                  beachmat_2.6.4               
##  [53] globals_0.15.0                goftest_1.2-3                
##  [55] tidygraph_1.2.1               rlang_1.0.2                  
##  [57] splines_4.0.4                 lazyeval_0.2.2               
##  [59] spatstat.geom_2.4-0           yaml_2.2.1                   
##  [61] reshape2_1.4.4                abind_1.4-5                  
##  [63] httpuv_1.6.5                  qvalue_2.22.0                
##  [65] tools_4.0.4                   ellipsis_0.3.2               
##  [67] spatstat.core_2.4-4           jquerylib_0.1.4              
##  [69] RColorBrewer_1.1-3            ggridges_0.5.3               
##  [71] Rcpp_1.0.8.3                  plyr_1.8.7                   
##  [73] sparseMatrixStats_1.2.1       progress_1.2.2               
##  [75] zlibbioc_1.36.0               RCurl_1.98-1.6               
##  [77] prettyunits_1.1.1             rpart_4.1-15                 
##  [79] deldir_1.0-6                  viridis_0.6.2                
##  [81] pbapply_1.5-0                 cowplot_1.1.1                
##  [83] zoo_1.8-10                    ggrepel_0.9.1                
##  [85] cluster_2.1.0                 magrittr_2.0.1               
##  [87] RSpectra_0.16-1               data.table_1.14.2            
##  [89] scattermore_0.8               DO.db_2.9                    
##  [91] lmtest_0.9-40                 RANN_2.6.1                   
##  [93] fitdistrplus_1.1-8            hms_1.1.1                    
##  [95] patchwork_1.1.1               mime_0.12                    
##  [97] evaluate_0.15                 xtable_1.8-4                 
##  [99] gridExtra_2.3                 compiler_4.0.4               
## [101] shadowtext_0.0.8              KernSmooth_2.23-18           
## [103] crayon_1.5.1                  htmltools_0.5.2              
## [105] ggfun_0.0.6                   mgcv_1.8-33                  
## [107] later_1.3.0                   tidyr_1.2.0                  
## [109] DBI_1.1.2                     ExperimentHub_1.16.1         
## [111] tweenr_1.0.2                  dbplyr_2.1.1                 
## [113] rappdirs_0.3.3                MASS_7.3-56                  
## [115] babelgene_22.3                Matrix_1.4-1                 
## [117] cli_3.2.0                     igraph_1.3.0                 
## [119] pkgconfig_2.0.3               rvcheck_0.2.1                
## [121] plotly_4.10.0                 spatstat.sparse_2.1-1        
## [123] bslib_0.3.1                   XVector_0.30.0               
## [125] yulab.utils_0.0.4             digest_0.6.27                
## [127] sctransform_0.3.3             RcppAnnoy_0.0.19             
## [129] spatstat.data_2.2-0           rmarkdown_2.14               
## [131] leiden_0.4.2                  fastmatch_1.1-3              
## [133] uwot_0.1.11                   DelayedMatrixStats_1.12.3    
## [135] curl_4.3                      shiny_1.7.1                  
## [137] lifecycle_1.0.1               nlme_3.1-152                 
## [139] jsonlite_1.8.0                BiocNeighbors_1.8.2          
## [141] limma_3.46.0                  viridisLite_0.4.0            
## [143] fansi_1.0.3                   pillar_1.7.0                 
## [145] lattice_0.20-41               fastmap_1.1.0                
## [147] httr_1.4.3                    survival_3.3-1               
## [149] GO.db_3.12.1                  interactiveDisplayBase_1.28.0
## [151] glue_1.4.2                    png_0.1-7                    
## [153] BiocVersion_3.12.0            bit_4.0.4                    
## [155] ggforce_0.3.3                 stringi_1.7.6                
## [157] sass_0.4.1                    blob_1.2.3                   
## [159] AnnotationHub_2.22.1          BiocSingular_1.6.0           
## [161] memoise_2.0.1                 irlba_2.3.5                  
## [163] future.apply_1.9.0